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Abstract 

The previously reported neBEM solver has been used to solve electrostatic problems 
having three-dimensional edges and corners in the physical domain. Both rectangular and 
triangular elements have been used to discrctize the geometries under study. In order to 
maintain very high level of precision, a library of C functions yielding exact values of poten- 
tial and flux influences due to uniform surface distribution of singularities on flat triangular 
and rectangular elements has been developed and used. Here we present the exact expres- 
sions proposed for computing the influence of uniform singularity distributions on triangular 
elements and illustrate their accuracy. We then consider several problems of electrostatics 
containing edges and singularities of various orders including plates and cubes, and L-shaped 
conductors. We have tried to show that using the approach proposed in the earlier paper on 
neBEM and its present enhanced (through the inclusion of triangular elements) form, it is 
possible to obtain accurate estimates of integral features such as the capacitance of a given 
conductor and detailed ones such as the charge density distribution at the edges / corners 
without taking resort to any new or special formulation. Results obtained using neBEM 
have been compared extensively with both existing analytical and numerical results. The 
comparisons illustrate the accuracy, flexibility and robustness of the new approach quite 
comprehensively. 

Keywords: Boundary element method, triangular element, singularity, electrostatics, po- 
tential, flux, capacitance, charge density, corner, edge. 

1 Introduction 

One of the elegant methods for solving the Laplace / Poisson equations (normally an integral 
expression of the inverse square law) is to set up the Boundary Integral Equations (BIE) which 
lead to the moderately popular Boundary Element Method (BEM). In the forward collocation 
version of the BEM, surfaces of a given geometry are replaced by a distribution of point singu- 
larities such as source / dipole of unknown strengths. The strengths of these singularities are 
obtained through the satisfaction of a given set of boundary conditions that can be Dirichlet, 
Neumann or of the Robin type. The numerical implementation requires considerable care ^ 
because it involves evaluation of singular (weak, strong and hyper) integrals. Some of the no- 
table two-dimensional (while all the devices are 3D by definition, useful insight is often obtained 
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by performing a 2D analysis) and three-dimensional approaches used to evaluate the singular 
integrals are discussed in [U [2] and O [H El E] and the references in these papers. It is well- 
understood that many of the difficulties in the available BEM solvers stem from the assumption 
of nodal concentration of singularities which leads to various mathematical difficulties and to 
the infamous numerical boundary layers [H El |33] when the source is placed very close to the 
field point ([2] and references [4-6] therein). While mathematical singularities (that occur when 
the source and field points coincide) have been shown to be artifacts, several techniques have 
been used to remove difficulties related to physical or geometrical singularities (that occur when 
boundaries are degenerate, i.e., geometrically singular, or due to a jump in boundary conditions) 
such as gaussian quadrature integration, mapping techniques for regularization, bicubic trans- 
formation, nonlinear transformation and dual BEM techniques [8]. The last technique seems to 
be a popular one and capable of dealing with a relatively wide range of similar problems. 

Departing from the approaches mentioned in the above references and many more to be 
mentioned below, we had shown in an earlier paper \XQ\ that many of these problems can be 
eliminated or reduced if we adopt a new paradigm in which the elements are endowed with 
singularities distributed on them, rather than assuming the singularities to be concentrated at 
specific nodal points. Despite a large body of literature, closed form analytic expressions for 
computing the effects of distributed singularities are rare [111 I12j . complicated to implement 
and, often, valid only for special cases [13l[lll|15]. For example, in [11], the integration of the 
Green function to compute the influence of a constant source distribution is modified to an "n- 
plane" integration. The evaluation of this integration involves co-ordinate transformations and 
the resulting expressions are rather complicated. In [12J, the Gauss-Bonnet concept is used in 
which the panel is projected onto a unit sphere and the solid angle is determined from the sum 
of the induced angles. The procedure and the resulting expressions are neither simple, nor easy 
to implement in a computer code. In fact, possibly due to these difficulties, these approaches 
have remained relatively unpopular and even in very recent papers it is maintained that for 
evaluating the influence due to source distributed on triangular elements in a general case, 
one must apply non-analytic procedures [14j . Thus, for solving realistic but difficult problems 
involving, for example, sharp edges and corners or thin or closely spaced elements, introduction 
of special formulations (usually involving fairly complicated mathematics, once again) becomes 
a necessity [HI l33l I16j . These drawbacks are some of the major reasons behind the relative 
unpopularity of the BEM despite its significant advantages over domain approaches such as 
the finite-difference and finite-element methods (EDM and EEM) while solving non-dissipative 
problems [IZKIl]. 

The Inverse Square Law Exact Solutions (ISLES) library developed in conjunction with the 
nearly exact BEM (neBEM) solver [lOj, in contrast, is capable of truly modeling the effect of 
distributed singularities precisely and, thus, is not limited by the proximity of other singular 
surfaces or their curvature or their size and aspect ratio. The library consists of analytic solutions 
for both potential and flux due to uniform distribution of singularity on flat rectangular and 
triangular elements. These close-form exact solutions, termed as foundation expressions, are in 
the form of algebraic expressions that are long but without complications and are fairly straight- 
forward to implement in a computer program. In deriving these foundation expressions, while 
the rectangular elements were allowed to be of any arbitrary size |10tll9|. the triangular element 
was restricted to be a right angled triangle of arbitrary size j20tl21t[22]. Since any real geometry 
can be represented through elements of the above two types (or by the triangular type alone), this 
library has allowed us to develop the neBEM solver that is capable of solving three-dimensional 
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potential problems involving arbitrary geometry. It may be noted here that any non-right- 
angled triangle can be easily decomposed into two right-angled triangles. Thus, the right-angled 
triangles considered here, in fact, can take care of any three-dimensional geometry. 

A set of particularly difficult problems to be dealt with by BEM is one that contains corners 
and edges and, in this work, we will attempt to solve several problems belonging to this set. The 
perfectly conducting bodies studied here are unit square plate, L-shaped plate, cube, L-shaped 
volume and two rectangular plates meeting at various angles and creating an edge. Besides being 
interesting and difficult, these solutions can have significant applications in micro electromechan- 
ical systems (MEMS), nano-devices, atomic force microscopy (AFM), electro-optical elements, 
micro-pattern gas detectors (MPGD) and many other disciplines in science and technology. For 
these problems, it is important to study integral features such as the capacitance of the conduc- 
tors, as well as detailed features such as the charge density, potential and flux on various surfaces 
of these objects including regions close to the geometric singularities. While several approaches 
including finite-difference method (EDM), finite element method (EEM), BEM and its variants 
such as the surface charge method (SCM) and various implementations of the Monte-Carlo 
technique (often coupled with Kelvin transformation) have been used to study these problems, 
only the latter two approaches, namely, BEM and Monte-Carlo technique are found to possess 
the precision necessary to model the curiously difficult electrostatics with acceptable levels of 
accuracy [23] . The volume discretization methods are known to be unsuitable because of the 
open nature of the problem and the inadequate representation of edge and corner singularities. 
Methods using Kelvin inversion (or quadratic inversion), although accurate, have been found 
incapable of handling planar problems. It may be noted here that despite the usefulness of 2D 
analysis, there are an overwhelming number of problems that need to be addressed in 3D. As 
a result, several interesting approaches have been developed to analyze edge and corner related 
problems in complete three dimensions, without even the assumption of axial symmetry. In 
order to maintain applicability in the most general scenarios, in this work we will deal with the 
problems of edge and corner as truly three-dimensional objects even when comparing the results 
with 2D analytic ones. 

The problem of estimation of capacitance of square plate and cube raised to unit volt has 
been studied by an especially large number of workers using entirely different approaches. In 
fact, these have been considered to be some of the major unsolved problems of electrostatics, 
of which a solution is said to have been given by Dirichlet and subsequently lost. One of 
the more popular numerical approaches used to explore these problems is the BEM / SCM 
[24:\ [25l [26l [T3l [271 [28] . Some of the solution attempts are more than a century old and yielded 
quite acceptable results. The later studies [131 [271 EH]) used the mesh refinement technique 
and extrapolation of N (the number of elements used to discretize a given body) to infinity 
in order to arrive at more precise estimates of the capacitance. In order to carry out this 
extrapolation, uniform charge density scenario has been maintained so that the form of charge 
distribution on individual segments becomes irrelevant. According to the [28], it is justified to 
use uniform charge density on individual elements because increase in complexity through the 
use of non-uniform charge density ultimately does not lead to computational advantage. In [28] . 
the author mentions that for the cube, the element sizes are made such that the charge on each 
element remains approximately a constant (independent of its distance from an edge) since this 
arrangement is found to give the most accurate results. 

The problem of estimation of the order of singularities at edges and corners of different 
nature is strongly coupled with the problem of estimation of integral properties such as the 
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capacitance of conductors of various shapes. Thus, this problem, which also has importance in 
relation to other areas of science and technology as discussed earlier, has attracted the attention 
of a large number of workers as well. Here, fortunately, some analysis has been possible using 
purely theoretical tools [MIEO], at least for two-dimensional cases. In [3T], the authors used a 
singular perturbation technique to obtain the singularity index at inside and outside corners of 
a sectorial conducting plate. Similarly, corner singularity exponents were numerically obtained 
in [28]. According to [32], it was possible to achieve accuracy of one in million through the use 
of FEM approximations for both electrodes and surface charge density, in addition to proper 
handling of edge and corner singularities. In this investigation, the Fichera's theorem was used 
to correctly describe the peculiarities of surface charge density behaviour in the vicinity of the 
electrode ribs and tips. 

According to another recent work [33j, low-order polynomials used to represent the corners 
and edges lead to errors in estimation of the derivatives of the potential, and that is the crux 
of the problem. Despite its advantages (no prior knowledge of singular elements and the order 
of singularity is necessary) and good convergence characteristics, the mesh refinement approach 
has been mentioned to be less accurate than two other methods, namely, 1) singular elements, 
and 2) singular functions. Among these, beforehand knowledge of the location and behaviour 
of the singularity in terms of the order of singularity is required for (1). The singular function 
approach (2) performs the best since it uses both the above and also information on singularity 
profile (corresponding to the eigenvalues and eigenvectors of the given geometry) . Unfortunately, 
this approach requires complex formulation and is generally restricted to 2D. Since the singular 
element approach uses the location information and only the order of the singularity, it is more 
versatile and popular. Thus, [33] implements the method of singular elements in 3D while 
attempting to use proper shape function / interpolation function to model correct singular 
behaviour at corners and edges. It is shown that singular elements are needed to accurately 
capture the behavior at singular regions, such as sharp corners and edges, where standard 
elements fail to give an accurate result. Unfortunately, singular elements must be defined before 
it is possible to apply either the singular function or the singular element approach. This is 
a nontrivial task since for a realistic device, there can easily be thousands, if not millions, of 
elements involved. According to [34], the manual classification of boundary elements based on 
their singularity conditions is an immensely laborious task if not outright impractical. In [Mj, the 
authors developed an algorithm to extract the regions where singularity arises by querying the 
geometric model for convex edges based on geometric information of the model. The associated 
nodes of the boundary elements on these edges were then retrieved and categorized according 
to different types of singularity configuration. The algorithm developed was implemented in the 
PATRAN command language (PCL) on the MSC/PATRAN platform [35] which is an industry 
standard finite-element pre- and post-processor allowing a high degree of customization. In order 
to determine the order of singularity, two-dimensional results have been used directly [29^ [36] 
for edges under the assumption that the point lies sufficiently far from any corner. In fact, in 
|34) . any point other than the vertices has been considered to be far enough. This approach 
is quite unlikely to be very accurate, especially when the dimensions of the devices are small 
enough so that no point may even be considered to be far enough. Later, while discussing the 
results obtained using our proposed approach, we will return to this issue again. A study on the 
effect of bias ratio (a ratio of the largest element length to the smallest element length, similar 
to r- mesh refinement) has led the authors [33j to state that while for normal elements this ratio 
should be around 4:1 for a given problem, the singular element approach works better with lower 
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bias ratio. 

Besides the above procedures to identify the singular elements and their orders, attempts 
have also been made to place nodes at optimal positions while solving problems involving edges 
and corners. For example, in fracture mechanics [37j, for handling singularities of the order 
0.5, mid-side nodes in quadratic elements is normally shifted to quarter point positions. Several 
other techniques such as least square, constrained displacement of side nodes adjacent to a crack- 
tip [38] have also been used to find the optimum position of the mid-side nodes for handling 
singularities of other orders. In [39], the authors have subsequently shown that this approach of 
shifting mid-side nodes is likely to produce singularities of order 0.5 only and it is not suitable to 
impose arbitrary singularity in isoparametric elements by simply shifting side nodes to assumed 
positions. 

In addition to BEM approaches, work has also been carried out using radically different 
approaches, for example, the Monte Carlo methods. An impressive array of work exist providing 
very accurate estimates of capacitance and variation of charge density near edges and corners 
[ini SU m [13 [13 [23 US]- According to [23], the use of BEM introduces unnatural estimate of 
the charge density distribution - not for shapes with smooth contours (discs or spheres), but for 
plates and cubes. It has been mentioned that the situation is worst for the corner singularities 
of plates in which case there are no other surfaces present (as in a cube) to weaken the order of 
singularity at the corner. We will discuss this issue when results for square plates are presented 
below. The other point is that the BEM can not satisfy the boundary conditions at or, at least, 
close to the edge because, the collocation points in BEM does not match the boundary of the 
device. In fact, it is mentioned that the solutions become unstable when the collocation points 
are shifted away from the centroid of the elements. This notion, as mentioned in the earlier 
paragraph, is not without counter-examples. Moreover, this is a point that we plan to take 
up in a future communication where we hope to illustrate that for the proposed formulation, 
shifting the collocation points to non-centroid locations does not lead to numerical instabilities. 
In addition, the approach of extrapolating capacitance has been criticized on the ground that 
they do not match for different amounts of shift. It has been mentioned that no formal error 
analysis exists for methods other than EDM and the extrapolation is purely empirical in nature. 
It has been observed that the apparent high accuracy may be illusory in nature and citing [47j . 
it has been emphasized that situation can become even worse by attempting extrapolation of 
results obtained using non-equivalent meshes. 

By developing a model that incorporates the truly distributed nature of sources / doublets / 
vortices on surfaces of three-dimensional geometries, we have recently shown jlO ^ 119 1 1^2] that it 
is possible to use the same formulation for studying a very wide range of problems (multi-scale, 
involving multiple layers of dielectric materials) governed by the Poisson's equation. Recently, 
we have extended the new formulation with the capability of including triangular elements as 
an option for discretizing arbitrary three-dimensional bodies [20^ I21j . Here, we present the 
expressions to evaluate the exact values of potential and fluxes at any arbitrary point due to 
uniform singularity distributed on right-angled triangular element. These expressions have been 
included as additional functions of the ISLES library and subsequently used in the neBEM 
solver in the manner usual to, probably, most of the BEM solvers available. Besides presenting 
the expressions, we have also shown results to illustrate the accuracy of the expressions under 
various circumstances. 

In this report, we have presented studies on the electrostatic configuration of several three- 
dimensional bodies, all of which contain corners and / or edges. The classic benchmark problems 
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of estimating the capacitance of a unit square plate and unit cube raised to unit volt have been 
addressed. Electrostatics of generic shapes such as L-shaped plate and L-shaped 3D conductors 
has also been analyzed. In all the above cases, the order of the singularity distribution near the 
edges and corners has been estimated in addition to estimating the capacitance of the conductors 
themselves. The singularity distributions obtained have been compared with theoretical and 
numerical studies carried out by earlier workers. The variation of the singularity along an edge 
between two corners at its ends has been studied, probably for the first time. Finally, the well- 
known problem of electrostatic configuration of two planes intersecting at different angles has 
been addressed. The two-dimensional counter-part of this problem is known to have analytic 
solution and has even been discussed in several textbooks on electromagnetics [291 [30]. Although 
accessible analytically, this problem seems to have been rarely solved using numerical techniques 
[48]. This benchmark problem is known to be a difficult one and, in order to test the proposed 
approach under difficult circumstances, we have computed electrostatic properties of a three- 
dimensional analogue of this problem close to the point of intersection for a very wide range 
of angle of intersection. Following the above studies, we have come to the conclusion that the 
proposed approach is capable of solving critical multi-scale problems governed by the Poisson's 
equation in a rather straight-forward manner. While higher bit accuracy, improved evaluation 
of transcendental functions, adaptive mesh generation and parallelization is expected to be of 
significant help, no special mathematical treatment or new formulation has been found to be 
necessary to deal with problems involving corners / edges and extremely closely spaced surfaces. 

According to ^S], using BEM, it is difficult to obtain physically consistent results close to 
these geometric singularities. Wild variations in the magnitude of the charge density have been 
observed with the change in the degree of discretization, the reason once again being associated 
with the nodal model of singularities. In contrast, using neBEM, we have obtained very smooth 
variation close to corner. Presence of oscillations seemingly acceptable to [23], has not occurred. 
In fact, oscillations close to edges and corners considered in this work seem to indicate numerical 
inaccuracy and have been treated accordingly. In addition to the shape, the magnitudes of the 
charge density have been found to be consistently converging to physically realistic values. These 
results clearly indicate that since the foundation expressions of the solver are exact, it is possible 
to find the potential and flux accurately in the complete physical domain, including the critical 
near-field domain using neBEM. In addition, since the singularities are no longer assumed to be 
nodal and we have the exact expressions for potential and flux throughout the physical domain, 
the boundary conditions no longer need to be satisfied at special collocation points such as 
the centroid of an element. Although consequences of this considerable advantage is still under 
study, it is expected that this feature will allow neBEM to yield even more accurate estimates for 
problems involving corners and edges since it should be possible to generate an over-determined 
system of equations by placing extra collocation points near the edges / corners. This will also 
allow the method to satisfy the boundary conditions of a given geometry at its true boundaries. 

It should be noted here that the exact expressions for triangular elements consist of a signif- 
icantly larger number of mathematical operations than those for rectangular elements presented 
in LlQJ- Thus, for any solver based on the ISLES library, it is more economical to use a mixed 
mesh of rectangular and triangular elements using rectangular elements as much as possible. 
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2 Governing equations and Exact Solutions 



In the following discussions, we will concentrate on the electrostatics of conducting geometries 
governed by the Poisson's equation. Using BEM approach, the Poisson's equation for electro- 
static potential 

V20(f) = -p{r)/eo 

can be solved to obtain the distribution of charges which leads to a given potential configuration. 
For a point charge q at r ' in 3D space, the potential (j){r) at r is known to be 



47reo|r — r'\ 

For a general charge distribution with charge density p{r'), superposition holds and results in 



47reo|r — r'l 
where 

G{r, P) = ^- — ^ 

47reo|r — r'l 

is the free space Green's function for the Laplace operator in 3D with eo, the permittivity of 
free space. Similarly, the field for a general charge distribution can be written as 

E{r) = -V0 

leading to 

E{r) = -V (^J G{r,r')p{r')dv' 
and, finally to, 

E{f) = f Pi^^'^i^-^y^' (2) 

The charge distribution can be obtained from equation ([1]) or ([2]) by satisfying the boundary 
conditions at collocation points known either in the form of potential (Dirichlet) or flux (Neu- 
mann) or a mixture of these two (Mixed/Robin) on material boundaries/surfaces present in the 
domain. 

Considering the Dirichlet problem only at present (for ease of discussion), the following 
integral equation of the first kind can be set up. 

G(f, r')p{r')dv' (3) 

vol 

In the above equation, (/)(r) is the potential at a point r in space and p{r') is the charge density 
at an infinitesimally small volume dv' placed at r' . The problem is, generally, to find p{r') as a 
function of space resulting the known distribution of 0(r). Once the charge distribution on the 
boundaries and all the surfaces is known, potential and field at any point in the computational 
domain can be obtained using the same equation ([3]) and its derivative. 

The primary step of the BEM technique is to discretize the boundaries and surfaces of a 
given problem. The elements resulting out of the discretization process are normally rectangular 
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or triangular though elements of other shapes are also used. Elements of triangular shape can 
be used to model geometries of any variety and, thus, is one of the most commonly used in many 
approaches of numerical simulation including FEM and BEM. In the collocation approach, the 
next step is to find out charge distribution on the elements that satisfies equation ^ following the 
given boundary conditions. The charge distribution is normally represented in terms of known 
basis functions with unknown coefficients. For example, in zero-th order formulations using 
constant basis function, which is also the most popular one among all the BEM formulations 
because of a good optimization between accuracy and computational complexity, the charge 
distribution on each element is assumed to be uniform and equivalent to a point charge located 
at the centroid of the element. This is the method that is referred to as the usual BEM in the 
rest of the paper. However, diverse varieties of basis function have been exercised to develop 
many more BEM formulations in order to represent the charge distribution on an element more 
efficiently so as to enhance the accuracy of the method. Since the potentials on the surface 
elements are known from the given potential configuration, equation ([3]) can be used to generate 
algebraic expressions relating unknown charge densities and potentials at the centroid of the 
elements. One unique equation can be obtained for each centroid considering influences of all 
other elements including self influence and, thus, the same number of equations can be generated 
as there are unknowns. In matrix form, the resulting system of simultaneous linear algebraic 
set of equations can be written as follows 

K-p = (^ (4) 

where K is the matrix consisting of influences among the elements due to unit charge density 
on each of them, p represents a column vector of unknown charge densities at centroids of the 
elements and (j) represents known values of potentials at the centroids of these elements. Each 
element of this influence coefficient or capacity coefficient matrix, K is a direct evaluation of 
an equation similar to equations ([I]) or ([2]) which represents the effect of a single element on a 
boundary/surface (obtained through discretization) on a point where a boundary condition of the 
given problem is known. While, in general, this should necessitate an integration of the Green's 
function over the area of the element, this integration is avoided in most of the BEM solvers 
through the assumption of nodal concentration of singularities with known basis function. The 
construction of the matrix implies that its diagonal elements are dominant through the presence 
of the Coulomb-type singularity in the kernel. This singularity has been shown to make the 
solutions well-defined in the class of rather smooth functions [32] . Since the right hand side of 
^ is known, in principle, it is possible to solve the system of algebraic equations and obtain 
surface charge density on each of the element used to describe the conducting surfaces of the 
detector following 

Once the charge density distribution is obtained, equations ([T]) and ([2]) can be used to obtain 
both potential and field at any point in the computational domain. 

Despite the elegance of formulation, the usual BEM suffers from several drawbacks that have 
resulted in its relative lack of popularity. Two of the most important ones can be mentioned 
as follows, (i) It is assumed that a surface distribution of charge density on an element can 
be represented by a nodal arrangement based on a chosen basis function, (ii) It is assumed 
that the satisfaction of the boundary condition at a predetermined point (or, through the use of 
known shape functions) is equivalent to satisfying the same on the whole element in a distributed 
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manner. The former assumption leads to infamous numerical boundary layer [8l [9l [33| [2] due to 
which the near-field solution in regions close to an element becomes erroneous. Thus the estima- 
tion of potential and field in near-field region close to the boundaries and surfaces by usual BEM 
is found to be inaccurate. This also leads to complications in solving problems involving closely 
spaced surfaces such as degenerate surfaces, edges, corners and other geometrical singularities. 
The degenerate surface refers to a boundary, two portions of which approach each other such 
that the exterior region between the two portions becomes infinitely thin. It is well known that 
the coincidence of two boundaries gives rise to an ill-conditioned problem. A number of special 
formulations has been developed to cope up with these problems but, unfortunately, most of 
these formulations are effective in a rather small subset of problems related to potential and 
field that are usually faced in reality. 

This problem has been resolved to a great extent through the development of the neBEM 
solver that uses exact integration of the Green's function and its derivative in its formulation. 
These integrations for rectangular and triangular elements having uniform charge density have 
been obtained as closed-form analytic expressions using symbolic mathematics [lUj . Thus they 
account for truly distributed nature of charge density on a given element. Besides the fun- 
damental change in the way the influence coefficient matrix is computed and the foundation 
expressions used for evaluating potential and field at any point after the charge density vector 
is solved for, most of the other features of neBEM are similar to any other BEM solver. 

The expressions for potential and flux at a point {X, Y, Z) in free space due to uniform source 
distributed on a rectangular flat surface having corners situated at (xi, zi) and (x2, Z2) has been 
presented, validated and used in [ini IS] and, thus, is not being repeated here. 

Here, we present the exact expressions necessary to compute the potential and flux due to a 
right-angled triangular element of arbitrary size, as shown in Fig(TJ It may be noted here that 
the length in the X direction has been normalized, while that in the Z direction has been allowed 
to be of any arbitrary magnitude, zm- From the figure, it is easy to see that in order to find 
out the influence due to triangular element, we have imposed another restriction, namely, the 
necessity that the X and Z axes coincide with the perpendicular sides of the right-angled triangle. 
Both these restrictions are trivial and can be taken care of by carrying out suitable scaling and 
appropriate vector transformations. It may be noted here that closed- form expressions for the 
influence of rectangular and triangular elements having uniform singularity distributions have 
been previously presented in [12\ I13j . However, in these works, the expressions presented are 
quite complicated and difficult to implement. In [lOj and in the present work, the expressions we 
have presented are lengthy, but completely straight-forward. As a result, the implementation 
issues of the present expressions, in terms of the development of the ISLES library and the 
neBEM solver, are managed quite easily. It is easy to show that the influence (potential) at 
a point P{X, Y, Z) due to uniform source distributed on a right-angled triangular element as 
depicted in FigH] can be represented as a multiple of 



in which we have assumed that xi = 0, zi = 0, X2 = 1 and Z2 = zm, as shown in the geometry 
of the triangular element. The closed-form expression for the potential has been obtained using 
symbolic integration [l9] which was subsequently simplified through substantial effort. It is found 
to be significantly more complicated in comparison to the expression for rectangular elements 
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Figure 1: Right-angled triangular element with x- length 1 and an arbitrary z-length, zm] P is 
the point where the influence (potential and flux) is being computed 

presented in p!0] and can be written as 
$ = 



Dii = ^{X-xiY + Y^ + {Z-z^)^- Z)i2 = ^{X-x^f + Y^ + {Z-Z2f 

D21 = ^J{X-X2f + Y^ + {Z-ZlY■, h = {X- xi) \Y\ ; h = {X - X2) \Y\ 
Si = sign{zi - Z); Ri=Y^ + {Z - zif 
Ei = {X + ZM^ - zmZ); E2 = iX-l- zmZ), 
G = zm{X + Z; Hi = Y^ + G{Z - zm); H2 = Y^ + GZ 




where, 



LPi 



G-izM\Y\ 



1 



log{ 



(ffi +GDi2)+i\Y\{Ei-izMDi2) 
-X + i\Y\ 
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^"^^ = G + ^ZM\Y\ ^''^ ^X^\ 

LP2 = 



1 , AH2 +GD2i)+i\Y\{E2-izMD2i) . 

Izm\Y\ l-X + i\Y\ 



G — izM 

1 , AH2 +GD2i)-i\Y\{E2-izMD2i), 

= gTJ^I^'^^ T^x^l ) 

and C denotes a constant of integration. 

Similarly, the flux components due to the above singularity distribution can also be repre- 
sented through closed-form expressions as shown below: 

ox 



dx 

^ ( (G)(LPi + LMi - LP2 - LM2) - i \Y\ {zm){LPi - LMi - LP2 + LM2) 
.^=^log(^^fl^)) + C (7) 

VI + ZM'^ V 1 + ZM'^D2l - E2 



-1 



{2zmY){LPi + LMi - LP2 - LM2) + i \Y\ {Sn{Y)G){LPi - LA'h - LP2 + LM2) 
+iS,Sn{Y){tanh-\?^^) - tanh-\^^) - tanh-\?^^) + tanh-\^-^)) ) + C 

(8) 

and, 

where Sn{Y) implies the sign of the Y-coordinate and C indicates constants of integrations. 
It is to be noted that the constants of different integrations are not the same. In addition 
to being extremely useful in the mathematical modeling of physical processes governed by the 
inverse square laws, these expression are expected to be useful as benchmark expressions for 
other approximate formulations. Being exact and valid throughout the physical domain, they 
can be used to formulate versatile solvers to solve multi-scale multi-physics problems governed 
by the Laplace / Poisson equations involving Dirichlet, Neumann or Robin boundary conditions. 
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3 Results and Discussions 



3.1 Exact expressions 



The expressions for the rectangular element have been validated in detail in [10]. Here, we present 
the results for triangular elements in fair detail. In FigOl we have presented a comparison of 
potentials evaluated for a unit triangular element by using the exact expressions, as well as 
by using numerical quadrature of high accuracy. The two results are found to compare very 
well throughout. Please note that contours have been obtained on the plane of the element, 
and thus, represents a rather critical situation. Similarly, Figl3] shows a comparison between 
the results obtained using closed-form expressions for flux and those obtained using numerical 
quadrature. The flux considered here is in the Y direction and is along a line beginning from 
(—2, —2, —2) and ending at (2, 2, 2). The comparison shows the commendable accuracy expected 
from closed form expressions. In Figs 4(a) and |4(b)] the surface plots of potential on the element 
plane (XZ plane) and y-flux on the XY plane have been presented from which the expected 
significant increase in potential and sharp change in the flux value on the element is observed. 
Thus, by using a small fraction of computational resources in comparison to those consumed in 
numerical quadratures, ISLES can compute the exact value of potential and flux for singularities 
distributed on triangular elements. 
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Figure 2: Potential contours on a triangular element computed using exact expressions and by 
numerical quadrature 
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Figure 3: Comparison of flux (in the Y direction) as computed by ISLES and numerical quadra- 
ture along a diagonal line 



3.1.1 Near-field performance 

In order to emphasize the accuracy of ISLES, we have considered the following severe situations 
in the near-field region in which it is observed that the quadratures can match the accuracy of 
ISLES only when a high degree of discretization is used. Please note that in these cases, the 
value of zm has been considered to be 10. In Figl^we have presented the variation of potential 
along a line on the element surface running parallel to the Z-axis of the triangular element (see 
Fig{T]) and going through the centroid of the element. It is observed that results obtained using 
even a 100 x 100 quadrature is quite unacceptable. In fact, by zooming on to the image, it can 
be found that only the maximum discretization yields results that match closely to the exact 
solution. It may be noted here that the potential is a relatively easier property to compute. 
The difficulty of achieving accurate flux estimates is illustrated in the two following figures. The 
variation of flux in the X-direction along the same line as used in Figl2]has been presented in 
Figini Similarly, variation of y-flux along a diagonal line (beginning at (-10,-10,-10) and ending 
at (10,10,10) and piercing the element at the centroid) has been presented in FiglTl From these 
figures we see that the flux values obtained using the quadrature are always inaccurate even if the 
discretization is as high as 100 x 100. We also observe that the estimates are locally inaccurate 
despite the use of very high amount of discretization (200 x 200 or 500 x 500). Specifically, in 
the latter figure, even the highest discretization can not match the exact values at the peak, 
while in the former only the highest one can correctly emulate the sharp change in the flux 
value. It is also heartening to note that the values from the quadrature using higher amount of 
discretization consistently converge towards the ISLES values. 
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(a) Potential surface 
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Figure 4: (a) Potential surface due to a triangular source distribution on the element plane, (b) 
Flux (in the Y direction) surface due to a triangular source distribution on the XY plane at Z=0 

3.1.2 Far field performance 

It is expected that beyond a certain distance, the effect of the singularity distribution can 
be considered to be the same as that of a centroidally concentrated singularity or a simple 
quadrature. The optimized amount of discretization to be used for the quadrature can be 
determined from a study of the speed of execution of each of the functions in the library and has 
been presented separately in a following sub-section. If we plan to replace the exact expressions 
by quadratures (in order to reduce the computational expenses, presumably) beyond a certain 
given distance, the quadrature should necessarily be efficient enough to justify the replacement. 
While standard but more elaborate algorithms similar to the fast multipole method (FMM) [50] 
along with the GMRES [CT] matrix solver can lead to further of computational efficiency, the 
simple approach as outlined above can help in reducing a fair amount of computational effort. 
In the following, we present the results of numerical experiments that help us in determining 
the far-field performance of the exact expressions and quadratures of various degrees that, in 
turn, help us in choosing the more efficient approach for a desired level of accuracy. 

In FiglH] we have presented potential values obtained using the exact approach, 100 x 100, 
10 X 10 and no discretization, i.e., the usual BEM approximation while using the zeroth order 
piecewise uniform charge density assumption. The potentials are computed along a diagonal line 
running from (-1000, -1000, -1000) to (1000, 1000, 1000) which pierces a triangular element of 
zm = 10. It can be seen that results obtained using the usual BEM approach yields inaccurate 
results as we move closer than distances of 10 units, while the 10 x 10 discretization yields 
acceptable results up to a distance of 1.0 unit. In order to visualize the errors incurred due 
to the use of quadratures, we have plotted FiglH where the errors incurred (normalized with 
respect to the exact value) have been plotted. From this figure we can conclude that for the 
given diagonal line, the error due to the usual BEM approximation falls below 1% if the distance 
is larger than 20 units while for the simple 10 x 10 discretization, it is 2 units. It may be mentioned 
here that along the axes the error turns out to be significantly more [lOj and the limits need 
to be effectively doubled to achieve the accuracy for all cases possible. Thus, for achieving 1% 
accuracy, the usual BEM is satisfactory only if the distance of the influenced point is five times 
the longer side of an element. Please note here that the error drops to 1 out of 10^ as the 
distance becomes fifty times the longer side. Besides proving that the exact expressions work 
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Figure 5: Variation of potential along a centroidal line on the XZ plane parallel to the Z axis 
for a triangular element: comparison among values obtained using the exact expressions and 
numerical quadratures 



equally well in the near- field as well as the far-field, this fact justifies the usual BEM approach 
for much of the computational domain leading to substantial savings in computational expenses. 



3.1.3 Comparison with multipole expressions 

We also compared the value of potential with estimates for the barycenter and other field point 
values as obtained from [14^ I15j. the expression having been slightly modified. At the barycen- 
ter, the expression given is exact, whereas, the multipole (monopole+quadrupole) expression is 
expected to be valid at any arbitrary location. The accuracy depends on the distance from the 
element as discussed in the papers, an approximate rule being that the accuracy is of the order 
of 10~^ when the distance of the field point from the barycenter of the triangle is more than 5.5 
times the longer side of the right triangle. For points less distant than the mentioned value, the 
triangle needs to be further segmented. Consulting Fig. [9] we may note that along a diagonally 
intersecting line, the usual BEM achieves this accuracy only when the distance of a field point 
is 8 times the larger side. The 10 x 10 discretization is as good for a field point that is distant 
twice the longer side. 

In Table [H we show the values estimated at the barycenter by ISLES, analytic [H] and 
numerical quadrature of different discretizations. Triangular elements of different sizes have 
been used keeping the x-side always of unit length. It is interesting to note that the right 
triangle for which two perpendicular sides are of equal length, it is most difficult to obtain the 
precise value of potential using quadrature. In fact, with a discretization of 2000 x 2000, we 
obtained the value of 2.407462, while with 5000 x 5000 we obtained 2.407323. 

In FigllO^ results from the above multipole expansion are compared with those obtained 
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Figure 6: Variation of flux in the X direction along a line on the XZ plane parallel to the Z axis 
for a triangular element: comparison among values obtained using the exact expressions and 
numerical quadratures 



Table 1: Comparison of estimated potential by ISLES, analytic and various quadratures 



ZM 


ISLES 


Analytic 


10 X 10 


100 X 100 


500 X 500 


0.1 


0.545069 


0.545069 


0.5410382 


0.5450810 


0.5450695 


1.0 


2.407320 


2.407320 


2.460945 


2.411947 


2.408161 


10.0 


5.450690 


5.450690 


5.257222 


5.450847 


5.450696 



using ISLES and numerical quadrature with different levels of discretization along a line passing 
from (-10,-10,-10) to (10,10,10), the range being reduced to a distance of -2 to +2 for ease 
of viewing. Similar comparison has been carried along a line parallel to the X-axis passing 
through the barycenter of the element and in the same plane as the element in Figllll As 
indicated correctly in |14^ I15j. the multipole expansion works fine at distances far enough. At 
close distances, only the ISLES results are acceptable. Other options exist in terms of further 
discretization. In that sense, both multipole expansion and numerical quadrature are likely to 
work but it may be difficult to decide on the required level of discretization a priori. Moreover, 
any advantage in terms of computational expenses may be lost due to the necessity of increased 
discretization. 



3.1.4 Speed of execution 

The time taken to compute the potential and fiux is an important parameter related to the overall 
computational efficiency of the codes. This is true despite the fact that, in a typical simulation, 
the time taken to solve the system of algebraic equation is far greater than the time taken to 
build the influence coefficient matrix and post-processing. Moreover, the amount of time taken 
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Figure 7: Comparison of flux (in the Y direction) along a diagonal line piercing the triangular 
element at the centroid: comparison among values obtained using the exact expressions and 
numerical quadratures 



to solve the system of equations tend to increase at a greater rate than the time taken to complete 
the other two. Thus, in fact, evaluation of longer expressions in the pre- and post-processing 
phase should hardly influence the efficiency of a BEM solver adversely. Moreover, due to the 
enhanced accuracy of the proposed foundation expressions, it will be possible to use significantly 
less number of elements to represent a given device, ultimately leading to much faster execution 
for a required level of accuracy. It should be mentioned here that the time taken in each of these 
steps can vary to a significant amount depending on the algorithm of the solver. In the present 
case, the system of equations has been solved using lower upper decomposition using the well 
known Grout's partial pivoting. Although this method is known to be very rugged and accurate, 
it is not efficient as far as number of arithmetic operations, and thus time, is concerned. It is also 
possible to reduce the time taken to pre-process (generation of mesh and creation of inffuence 
matrices), solve the system of algebraic equations and that for post-process (computation of 
potential, fiux at required locations, estimation of other properties such as capacitance, force) 
by adopting faster algorithms, including those involving parallelization. 

In order to optimize the time taken to generate the inffuence coefficient matrix and that to 
carry out the post-processing, we carried out a small numerical study to determine the amount 
of time taken to complete the various functions being used in ISLES, especially those being used 
to evaluate the exact expressions and those being used to carry out the quadratures. The results 
of the study (which was carried out using the linux system command gprof) has been presented 
in the following Table [2j 

Please note that the numbers presented in this table are representative and are likely to 
have statistical fiuctuations. However, despite the fiuctuations, it may be safely concluded that 
a quadrature having only 10 x 10 discretization is already consuming time that is comparable 
to that needed for exact evaluation. Thus, the exact expressions, despite their complexity. 
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Figure 8: Potential along a diagonal through the triangular element computed using exact, 
100 X 100, 10 X 10 and usual BEM approach 



Table 2: Time taken to evaluate exact expression of ISLES, Usual BEM and various quadratures 

Method ISLES Usual BEM 10 x 10 100 x 100 500 x 500 

Time for 0.6 fis 25 ns 2 400 fis 10 ms 

rectangular 
element 

Time for 0.8 ij.s 25 ns 2 fis 400 fis 10 ms 

two triangular 

elements 



are extremely efficient in the near-field which can be considered at least as large as 0.5 times 
the larger side of a triangular element (please refer to FigE]). In making this statement, we 
have assumed that the required accuracy for generating the influence coefficient matrix and 
subsequent potential and flux calculations is only 1%. This may not be acceptable at all under 
many practical circumstances, in which case the near-field would imply a larger volume. 

Some of the advantages of using the ISLES library, based on the presented foundation ex- 
pressions, are mentioned below: 

• For a given level of discretization, the estimates are more accurate, 

• Effective efficiency of the solver improves, as a result, 

• Large variation of length-scales, aspect ratios can be tackled, 

• Thinness of members or nearness of surfaces does not pose any problem, 

• Curvature has no detrimental effect on the solution. 
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Figure 9: Error along a diagonal through the triangular element computed using 100 x 100, 
10 X 10 and usual BEM approach 



• The boundary condition can be satisfied anywhere on the elements, i.e., points other than 
the centroidal points can be easily used, if necessary (for a corner problem, may be), 

• The same formulation, library and solver is expected to work in majority of physical 
situations. As a result, the necessity for specialized formulations of BEM and associated 
complications can be greatly minimized. 

3.2 Electrostatics of two- and three-dimensional bodies having corners and 
edges 

3.2.1 Square plate and Cube 

The capacitance value estimated by the present method has been compared with very accurate 
results available in the literature (using BEM and other methods). The results obtained using 
the neBEM solver is found to be among the most accurate ones available till date as shown in 
Table [3l Please note that we have neither invoked symmetry nor used extrapolation techniques 
to arrive at our result presented in the table. 

In Table HI we compare potentials at the center and along an edge of the unit cube as 
obtained using neBEM with those from [32] in which the authors use analytical techniques to 
determine the order of singularity at the singular regions of a cube. From this table, we find 
that it has been possible for [32j to maintain accuracy of 10~^ in 1 for the potential values 
at the cube center, 10~^ to 10~^ in 1 along an edge and 10^^ in 1 along edge but close to a 
vertex of the cube. For similar locations, results using neBEM indicate that an accuracy of 
10~^ is maintained at the cube center, 10~^ at the edge and 10"'^ on the edge but close to the 
vertex. Thus, the proposed approach has been able to achieve accuracy that is comparable to 
those achieved by [32] that uses the Fichera's theorem to ensure proper variation of singularities 
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Figure 10: Potential along a diagonal through the triangular element computed using ISLES, 
multipole expansion, 100 x 100, 100 x 10 and usual BEM approach 



near edges and vertices. Interestingly, at critical locations near the cube vertex, results from 
neBEM are found to be better than [32] by a significant amount. We have been able to maintain 
this accuracy of 10~^ as close as upto Ijim to the cube vertex which is extremely encouraging. 
Unfortunately, we have not been able to compare our results with other numerical results at 
distances less than a mm from the vertex. It is encouraging to note, however, that while for 
|32j . the error is 4.8 x 10^'^ when the evaluation point is 1mm away from the vertex, neBEM 
commits an error of a similar amount only when the point is lO/um away from the vertex. The 
error for neBEM becomes larger by a small amount (0.4 x 10~^) only when the evaluation point 
is as close as a micron to the vertex. 

Next, we consider the problem of determining charge density distribution at corners and 
edges of the above geometries. Problems of this nature are considered to be challenging for 
any numerical tool and especially so for the BEM approach. In Table [5l we have presented the 
estimates of the order of singularity at the vertex or the edge as done by methods as diverse as 
singular perturbation [31] , BEM [28] , last-passage and walk on spheres [lU US] , application of 
Fichera's theorem |32] . singular element approach [3H] and the presented approach. From the 
table, it is clear that there is good agreement among all the methods. As in [28] properties on 
the element next to a corner or edge has been ignored while carrying out the least-square fits. 
Points were included in the fit as long as the maximum mismatch between the fitted line and 
the computed value was less than typically 1%, which also allows us to include points as long as 
they fall closely on a straight line. Following this approach, we could use values of singularities 
even upto 0.15 (for a plate or cube of unit side) from the relevant edge or corner while fitting 
the lines (see FigJ131 as an example where we used both triangular and rectangular elements for 
discretizing a square plate and obtained a singularity index of 0.7057 and 0.7068, respectively). 
These can be compared to the facts that in [^, two points next to a corner was excluded, as 
well as all points at distances beyond 0.05 from a corner. However, the exact value of the index 
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Figure 11: Potential along a line parallel to X through the centroid of a triangular element 
computed using ISLES, multipole expansion, 100 x 100, 10 x 10 and usual BEM approach 



is, to a certain extent, dependent on the discretization and the details of the least square fitting 
procedure. Thus, it may not be very prudent to attach great significance to the obtained values 
except noting that they agree with each other and also agree with the theoretical estimates, 
wherever available. Thus, it may be difficult to accurately ascertain the singularity index at 
an arbitrary corner or edge. This difficulty can lead to problems for methods that depend on 
beforehand knowledge of the order of singularity. 

In the following study, we have presented estimates of charge density very close to the flat 
plate corner as obtained using neBEM. This has been carried out to investigate the objection 
raised against the BEM in [23] where the author states that severe oscillations in the charge 
densities are expected close to the corner and edge of plates because the BEM cannot correctly 
model the edge / corner of physical devices. Please note that for this study, the boundary 
conditions have been satisfied at the centroids of each element although the neBEM has the 
capability of satisfying boundary conditions at locations other than the centroid. In Fig ll2l 
charge densities very close to the corner of the flat plate estimated by neBEM using various 
amounts of discretization have been presented. It can be seen that each curve follows the same 
general trend, does not suffer from any oscillation and is found to be converging to a single 
curve. This is true despite the fact that there has been almost an order of magnitude variation 
in the element lengths. Thus, we can safely conclude that the estimates obtained using neBEM 
do not suffer from the numerical instabilities mentioned in |23| . 

In FigllSt we present a least-square fitted straight line matching the charge density as ob- 
tained using the highest discretization in this study. Results using both triangular and rectan- 
gular elements are presented and it is found that the slope of the fitted line is 0.7057 when the 
elements are triangular, whereas, it is 0.7068 when we use rectangular elements. Both the values 
compare very well with both old and recent estimates of the order of singularity as shown in 
Table [5l In Fig fT^ we have shown how the slope of the fitted line changes along the edge of a 
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Table 3: Comparison of capacitance values 



Reference 


Method 


Plate (pF) / 4 ireo 


Cube (pF) / 4 ^eo 


[21] and [25] 


SCM 


0.3638 




m 


SCM 


0.362 


0.6555 


[52] 


SCM 


0.367 




[13] 


Refined SCM 


0.3667894 ± 1.1 x 10"^ 


0.6606747 ± 5 x IQ-'^ 




and Extrapolation 






m 


Refined BEM 


0.3667874 ± 1 x 10"^ 


0.6606785 ± 6 x 10"^ 




and Extrapolation 






m 


Numerical Path 


0.36684 


0.66069 




Integration 






[23] 


Random Walk 


0.36 ±0.01 


0.6606 ± 1 X 10"^ 


m 


Random Walk 




0.6606780 ± 2.7 x IQ-^ 


[33] 


Singular element 




0.6606749 


[28] 


Refined BEM 


0.3667896 ± 8 x 10"^ 


0.6606767 ± 4 x IQ-^ 




and Extrapolation 






m 


Robin Hood 




0.6606786 ± 8 x IQ-^ 




and Extrapolation 






This work 


neBEM 


0.3667524 


0.6606746 



Table 4: Comparison of potential at the center and along an edge of a unit cube 



X 


Y 


Z 


Exact 


[32] 


Error in [32j 


neBEM 


Error in neBEM 











1 


0.999990 


-1.0 X 10^^ 


1.000001 


1.0 X lO-'* 


0.4 


0.5 


0.5 


1 


0.9996 


-4.0 X 10"'' 


0.9994362 


-5.638 X 10-^ 


0.45 


0.5 


0.5 


1 


0.99986 


-1.4 X 10"'' 


0.9995018 


-4.982 X 10-^ 


0.49 


0.5 


0.5 


1 


1.0013 


1.3 X 10-3 


0.9991151 


-8.849 X 10"^ 


0.499 


0.5 


0.5 


1 


1.0048 


4.8 X 10-3 


0.9987600 


-1.24 X 10-3 


0.4999 


0.5 


0.5 


1 






0.9974398 


-2.56 X 10-3 


0.49999 


0.5 


0.5 


1 






0.9951335 


-4.8 X 10-3 


0.499999 


0.5 


0.5 


1 






0.9945964 


-5.4 X 10-3 



square plate as we move away from a corner of a square plate. From the figure, it is apparent 
that the change in the singularity index along an edge of a square plate, can be quite significant 
and only when we are in reality close to the middle of the edge, the analytic value of 0.5 can 
be used with confidence for the order of singularity. This observation is significant especially 
for the singular element and singular function methods where prior knowledge of this parameter 
plays a crucial role in determining the accuracy of the solution. 

3.2.2 L-shaped plate and volume 

Despite being geometries of generic nature, these L-shaped geometries have received relatively 
less attention. Here, we have estimated the capacitance and orders of singularity at various 
important locations, e.g., inner and outer corners and compared these values with available 
numerical results. It should be mentioned here that in |33[ I34j. the L-shaped volume has been 
described very nicely in relation to the varying nature of singularities it contains. Capacitance 
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Table 5: Comparison of estimation of order of singularity 



Reference Method 



Plate corner Plate edge Cube vertex Cube edge 



m 


Analytic 




0.5 




0.333 


m 


Numerical shooting 


0.7034 








[53] 


BEM 






0.5468 




m 


Walk on spheres 


0.7034 




0.5381 / 












0.5458 




m 


Surface Charge 


0.704 




0.540 




m 


Surface Charge 






0.558 


0.333 


m 


Singular element 






0.5475 




m 


Random Walk 










[32] 


Fichera's theorem 


0.7015 




0.5454 




m 


Walk on planes 


0.7034 




0.5457 





This work neBEM 0.7068 0.4994 0.5539 0.332 



of the L-shaped volume conductor has turned out to be 112.1497pF. The estimate matches 
extremely well with the value of 112.15pF that has been used as a reference value in [33]. In 
Table [H] we present the numerical values of the order of singularities at inner corners of the 
L-shaped plate and volume conductor. While the former matches well with |31j . the latter is not 
quite close to the estimate in [33], the difference being close to 20%. Finally, in FigfTSj we show 



Table 6: Comparison of estimation of order of singularity 



Reference 


Method 


Plate inner 


L volume inner 






corner 


corner 


[31] 


Numerical shooting 


0.1854 




[33] 


Singular element 




0.1104 


This work 


neBEM 


0.1840 


0.0896 



how the magnitude of charge density changes on a L-shaped plate. The remarkable difference 
between external and internal corners in terms of charge density concentration is very clearly 
observed in this figure. 

3.2.3 Edge problem having analytical solutions in 2D 

While none of the above problems have analytic solutions, a closely related problem has well- 
known analytic solution in the two-dimensional case [29]. In fact, while discussing the earlier 
problems, other workers and we have quite often referred to this solution obtained in standard 
textbooks as an exercise in the method of separation of variables. 

We have considered a three-dimensional equivalent of the geometry as presented in [29] in 
which two conducting planes intersect each other at an angle /?. The planes are assumed to be 
held at a given potential. In order to specify the boundary conditions conveniently, a circular 
cylinder is also included that just encloses the two intersecting plane, has its center at the 
intersection point and is held at zero potential [IS]. The general solution in the polar coordinate 
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Figure 12: Corner charge density estimated by neBEM using various sizes of triangular elements 



system (p, </>) for the potential {^) close to the origin in this problem has been shown to be 

oo 

^>(p, (t)) = V + Y^am /o'""/^ sm(m7r0//3) (10) 

m=l 



where the coefficients am depend on the potential remote from the corner at p = and V 
represents the boundary condition for ^ for all p > when = and (j) = (3. In the present case 
where a circular cylinder just encloses the two plates, the problem of finding out am reduces to 
a basic fourier series problem with a well known solution 

4 

am = foroddm (11) 

mvr 

It may be noted here that the series in (|10|) involve positive powers of p'^^^ , and, thus, close to 
the origin (i.e., for small p), only the first term in the series will be important. The electric field 
components {Ep, E^) are 

oo 

Ep{p, </>) = E «™ p(™"/^)~'sm(m7r0//3) (12) 

7Tl=l 

OO 

E^{p, 0) = E a„ p(™-/^)-icos(m7r</.//3) (13) 

'"^ m=l 

The surface charge densities (a) at = and (j) = (3 are equal and are approximately 

<P) = ^^^-%P'^"'-' (14) 
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Figure 13: Variation of charge density with increasing distance from the corner of the unit square 
plate and a least-square fitted straight line: slope of the fitted line is 0.7068 



Thus, the field components and the charge density near p = both vary with distance as 
and this fact is expected to be reflected in a correct numerical solution as well. 

While the above theoretical solution is a two-dimensional one, we have used the neBEM 
solver to compute a three-dimensional version of the above problem. In order to reproduce 
the two-dimensional behavior at the mid-plane, we have made the axial length of the system 
sufficiently long, viz., 10 times the radius of the cylinder. The radius of the cylinder has been 
fixed at one meter, while the length of the intersecting flat plates has been made a mm shorter 
than the radius to avoid the absurd situation of having two values of the voltage at the same 
spatial point. 

The cylinder has been discretized uniformly in the angular and axial directions. The flat 
plates have also been uniformly discretized in the axial direction. In the radial direction, however, 
the flat plate elements have been made successively smaller towards the edges using a simple 
algebraic equation in order to take care of the fact that the surface charge density increases 
considerably near the edges. From Tables [71 [8] and [9l we can compare the accuracy of neBEM 
results with other analytical and numerical results. The two ends of the range of angles, 360 
and 90, represent particularly difficult situations. The former is difficult due to the very large 
concentration of charge density close to a sharp edge and the resulting large electric field. The 
latter is difficult due to the fact that to truly simulate a null point in a concave corner, extremely 
precise estimates are necessary to ensure cancellation of electric field. Throughout the range, the 
neBEM results are found to be very accurate except at the location that is just l//m away from 
the edge. Even at this location, for all the convex edges, the results are reasonable and surely 
comparable to the only other numerical result available. The neBEM estimates are unacceptable, 
however, at locations less than tens of microns away from the null corner of the concave corner. 
It may be noted however, that the problem is not at all inherent to the formulation and is clearly 
related to the size of elements used in the vicinity of the corner. Here, on our desktop PC with 
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Figure 14: Variation of the slope of the fitted hnes along an edge of a square plate 



2GB RAM, we could use a spatial resolution of around a micron close to the corner despite using 
a large profiling factor. And this was at the steep cost of having elements with extremely large 
aspect ratios (1 : 10^). We believe that these are the factors that have resulted in the inaccuracy 
of the presented results when the locations considered were a micron away from the edge. 





Table 7: 


Electric field close to a 360deg edge 




Distance 


Analytical 


ELECTRO 


Error (%) 


neBEM 


Error (%) 


0.8 


0.3954180 


0.3954213 


0.00059 


0.3950786 


-0.086 


0.1 


1.830153 


1.830155 


0.00010 


1.830110 


-0.002 


0.01 


6.303166 


6.303172 


0.000094 


6.305784 


-0.041 


0.001 


20.11157 


20.11122 


0.0018 


20.11963 


-0.04 


0.0001 


63.65561 


63.64274 


0.020 


63.64780 


-0.012 


0.00001 


201.3148 


200.88 


0.22 


200.5488 


-0.3 


0.000001 


636.6191 


621.25 


2.4 


621.6034 


-2.36 



To end this section, we present one contour plot of electric field for a convex corner. From 
the Fig.??, it is evident that the intensity of the field increases only very close to the convex 
edge. Please note that the dimensions in the figure are mentioned in microns. 

4 Conclusion 

An efficient and robust library, ISLES, for solving potential problems in a large variety of science 
and engineering problems has been presented. Exact closed-form expressions used to develop 
ISLES have been validated throughout the physical domain (including the critical near-field 
region) by comparing these results with results obtained using numerical quadrature of high 
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Figure 15: Charge density estimated by neBEM on a L-shaped conducting plate raised to unit 
volt 





Table 8: 


Electric field close to a 270deg edge 




Distance 


Analytical 


ELECTRO 


Error (%) 


neBEM 


Error (%) 


0.8 


0.5246997 


0.524710 


0.0019 


0.5241510 


-0.105 


0.1 


1.747623 


1.747621 


0.00014 


1.747953 


-0.018 


0.01 


3.931433 


3.931284 


0.0038 


3.933242 


-0.046 


0.001 


8.487415 


8.4854 


0.023 


8.491335 


-0.046 


0.0001 


18.28732 


18.202 


0.46 


18.29270 


-0.029 


0.00001 


39.39902 


35.80 


9.1 


39.30955 


-0.227 


0.000001 


84.88264 


57.10 


32.7 


80.74309 


-4.877 



accuracy and with those obtained using quadrupole expressions. Algorithmic aspects of this 
development have also been touched upon. The neBEM solver that uses foundation expressions 
being evaluated by ISLES, has been used to solve several corner and edge electrostatic problems. 
Several classic benchmark problems such as those related to unit square plate and unit cube 
have been solved to very high precision. Charge density values at critical geometric locations like 
corners have been found to be numerically stable and physically acceptable. Values of singularity 
indices at different corners and edges have been estimated and compared with other analytic 
and numerical estimates. The agreement among the different approaches has been found to 
be quite acceptable. It has been observed that the variation of this index along the edge of 
a conducting body is non-negligible implying caution necessary for methods that need prior 
knowledge of these indices to solve a given problem. Finally, using the same solver, a three- 
dimensional equivalent of an edge problem has been solved for which analytic solution exists in 
two-dimensions. Detail comparison with this problem and those stated earlier have led us to 
believe that the neBEM approach is a precise, flexible and robust solver that works over a very 
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Figure 16: Electric field distribution very close to a convex edge. 



wide range of problems. Several advantages over usual BEM solvers and other specialized BEM 
solvers have been briefly mentioned. Some of the criticisms leveled against the BEM approach 
in general have been addressed in this work, while we expect to solve some of the remaining in 
future communications. Work is also under way to make the code more efficient through the 
implementation of faster algorithms and parallelization. 
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Table 


9: Electric 


field close to 90deg and 225de 


g edges 






Distance 


0.8 


0.1 


0.01 


0.001 


0.0001 


0.00001 


0.000001 


90 degrees 


Analytic 


1.445221 


2.546224e-l 


2.546479e-2 


2.546479e-3 


2.546479e-4 


2.546479e-5 


2.546479e-6 


neBEM 


1.444098 


2.547710e-l 


2.548018e-2 


2.547723e-3 


2.778723e-4 


3.7601 17e-3 


1.518850 


Error (%) 


-0.001123 


0.058361 


0.060436 


0.048846 


9.120200 


usable(?) 


inaccurate 


225 degrees 


Analytic 


0.6266090 


1.574802 


2.556973 


4.055022 


6.426876 


10.18592 


16.14359 


neBEM 


0.6259308 


1.575165 


2.557946 


4.056545 


6.428439 


10.16933 


15.63912 


Error (%) 


-0.108233 


-0.023050 


-0.038052 


-0.037558 


-0.024313 


-0.162871 


-3.124893 
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